/*------------------------------------------------------------------------
 *  Copyright (C) 2008  Luis M. de la Cruz
 *
 *  This file is part of TUNA::RBF
 *
 *  TUNA::RBF is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  TUNA::RBF is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ------------------------------------------------------------------------*/

template<class T_sfield>
void DiffOperators::vorticity(T_sfield &vort1, 
			      T_sfield &vort2,
			      T_sfield &vort3, 
			      const T_sfield &u, 
			      const T_sfield &v,
			      const T_sfield &w,
			      double dx, double dy, double dz) {

    int bi = u.lbound(firstDim) + 1, ei = u.ubound(firstDim) - 1;
    int bj = u.lbound(secondDim) + 1, ej = u.ubound(secondDim) - 1;
    int bk = u.lbound(thirdDim) + 1, ek = u.ubound(thirdDim) - 1;

    double du_dy, du_dz, dv_dx, dv_dz, dw_dx, dw_dy;

    for(int i = bi; i <= ei; ++i) 
	for(int j = bj; j <= ej; ++j) 
	    for(int k = bk; k <= ek; ++k) 
	    {
		du_dy = ( u(i,j+1,k  ) - u(i,j-1,k  ) ) / (2 * dy);
		du_dz = ( u(i,j  ,k-1) - u(i,j  ,k+1) ) / (2 * dz);

		dv_dx = ( v(i+1,j,k  ) - v(i-1,j,k  ) ) / (2 * dx);
		dv_dz = ( v(i  ,j,k+1) - v(i  ,j,k-1) ) / (2 * dz);

		dw_dx = ( w(i+1,j  ,k) - w(i-1,j  ,k) ) / (2 * dx);
		dw_dy = ( w(i  ,j+1,k) - w(i  ,j-1,k) ) / (2 * dy);

		vort1(i,j,k) = dw_dy - dv_dz;
		vort2(i,j,k) = du_dz - dw_dx;
		vort3(i,j,k) = dv_dx - du_dy;
	    }
}

template<class T_sfield>
void DiffOperators::vorticity(T_sfield &vort1,
			      const T_sfield &u, 
			      const T_sfield &v,
			      double dx, double dy) {

    int bi = u.lbound(firstDim) + 1, ei = u.ubound(firstDim) - 1;
    int bj = u.lbound(secondDim) + 1, ej = u.ubound(secondDim) - 1;

    double du_dy, du_dz, dv_dx, dv_dz, dw_dx, dw_dy;

    for(int i = bi; i <= ei; ++i) 
	for(int j = bj; j <= ej; ++j) 
	    {
		du_dy = ( u(i  ,j+1) - u(i  ,j-1) ) / (2 * dy);
		dv_dx = ( v(i+1,j  ) - v(i-1,j  ) ) / (2 * dx);
		vort1(i,j) = dv_dx - du_dy;
	    }
}





